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ABSTRACT 

The  objective  of  the  present  note  is  to  provide  the  basic  concepts  of  numerical  modelling  of  cavitation. 
After  a  brief  description  of  the  cavitation  phenomena,  the  modelling  conceps  are  discussed.  A  specific 
attention  is  given  to  the  development  of  a  barotropic  model  and  its  integration  in  a  3D  unsteady  CFD 
code.  Influence  of  numerical  and  physical  parameters  is  detailed  on  two  two-dimensional  configurations 
(hydrofoil  and  venturi).  Examples  of  calculations  on  various  geometries  ( inducers  and  impellers )  are 
detailed  to  validate  robustness  of  the  computational  methodology.  Modifications  of  internal  flows  due  to 
development  of  cavitation  are  discussed  and  an  analysis  of  the  mechanisms  leading  to  the  head  drop  is 
presented.  Simulations  of  cavitation  instabilities  in  blade  cascades  are  presented  and  physical 
mechanisms  are  discussed. 
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Uref  Reference  velocity  (m/s) 

V  Absolute  velocity 

yparoi  distance  to  the  wall  (m) 


ux  vitesse  de  frottement :  UT 

,  „+  _  PUr  yparoi 

y  distance  to  the  wall :  y  — 

M 

a  Void  fraction 

c  Turbulent  dissipation 

y  Surface  tension  (J/m2) 

p  Viscosity 

pt  Turbulent  viscosity 

p  Density  (kg/m3) 

pi,  pv  Liquid  and  vapour  density  (kg/m3) 

.  .  Vref-Pvap 

a  Cavitation  parameter:  <J= — - - 

\puif 

co  Rotation  speed  (rad/s) 

Q  Cell  volume  (m3) 

1.0  INTRODUCTION 

In  the  frame  of  the  development  of  a  rocket  engine  turbopump,  cavitation  is  one  of  the  most  limitating 
hydraulic  constraints.  For  this  reason,  CNES  and  Snecma  have  been  working  in  collaboration  with  LEGI 
laboratory  and  Numeca  in  order  to  develop  a  CFD  code  being  able  to  predict  the  cavitating  flows 
developing  in  rocket  engine  turbopumps.  For  this  presentation,  cavitation  will  be  considered  assuming  no 
thermal  effects. 

1.1  Cavitation 

Cavitation  relates  to  the  liquids  submitted  to  pressure  lower  than  their  saturating  vapour  pressure.  The 
liquid  phase  is  then  metastable  and  operates  a  transition  to  vapour  phase  by  formation  of  bubbles:  this 
process  is  called  cavitation.  When  cavitation  appears  in  a  flow,  the  physical  properties  are  deeply  modified 
in  particular  because  of  a  brutal  increase  in  compressibility  of  the  gas/liquid  mixture. 

1.1.1  Vapour  Structures 

From  a  macroscopic  point  of  view,  cavitation  can  appear  in  various  forms:  by  attached  cavities,  by 
bubbles  or  by  localised  vapour  structures  in  shear  layers.  The  internal  structure  of  the  cavities  is  complex 
and  still  badly  known.  However,  experimental  studies  such  as  those  performed  by  Stutz  [20]  have 
provided  major  informations  on  the  characterisation  of  void  fraction  and  velocity  distributions.  Thus,  it 
has  been  observed  that  the  cavitation  sheets  were  very  heterogeneous  with  void  fraction  values  ranging 
from  more  than  80  %  in  the  upstream  part  of  a  stable  cavity  to  less  than  10%  in  its  downstream  part.  In  the 
major  part  of  the  unsteady  cavities  the  maximum  values  of  void  fraction  do  not  exceed  20%. 

1.1.2  Cavitation  and  Turbomachinery 

The  operation  of  an  hydraulic  machine  generates  localised  depressions  on  the  blades.  When  the  decrease 
of  static  pressure  is  sufficient,  it  can  lead  to  the  appearance  of  cavitation  that  may  deteriorate  the 
performances  of  the  machine.  Consequently,  cavitation  has  to  be  considered  as  an  operational  limit  of  the 
hydraulic  machines: 


3-2 


RTO-EN-AVT-1 43 


Numerical  Modelling  of  Cavitation 


•  in  term  of  maximum  mass  flow 

•  in  term  of  height  of  installing.  This  last  criterion  is  extremely  important  since  it  is  mainly 
responsible  for  the  realization  costs  in  projects  like:  hydroelectric  power  plant,  pressurisation 
devices  installed  on  the  water  aduction  networks,  pumping  stations. . . 


Figure  1.1 :  Modification  of  a  Pump  Characteristic  Curve. 

Finally,  we  focus  on  the  problems  that  have  motivated  our  work:  cavitation  developing  in  the  rocket 
engine  turbopumps. 


◄ - - ► 

Pump  Turbine 

Figure  1.2:  Vulcain  LH2  Turbopump. 


When  a  turbopump  is  running,  cavitation  may  occur  in  the  axial  stage  that  is  called  inducer.  Vapour 
cavities,  when  they  are  highly  developed,  cause  a  decrease  of  the  pressure  rise  provided  by  the  pump.  The 
head  drop  occurs  in  a  more  or  less  brutal  way  according  to  the  flow  rate.  In  general,  one  observes  that  the 
lower  the  mass  flow,  the  more  brutal  the  head  drop. 

The  second  consequence  of  cavitation  development  in  an  inducer  is  the  appearance  of  instabilities. 
Cavitation  plays  a  role  at  system  level  (modification  of  hydraulic  frequencies,  damping)  and  at  pump 
level.  We  focus  on  the  instabilities  that  are  intrinsic  to  the  pump.  Instabilities  are  associated  with  a  rupture 
of  the  flow  symmetry;  thus,  several  configurations  can  be  distinguished  according  to  the  cavity  lengths  on 
each  blade: 
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Figure  1.3:  Cavitation  Instabilities  [4]. 

Two  major  instability  configurations  have  been  observed: 

•  Alternate  cavitation  is  characterized  by  the  coexistence  of  two  small  cavities  intercalated  between 
two  longer  cavities. 

•  Super  synchronous  cavitation  is  characterized  by  the  presence  of  a  small  cavity  surrounded  by 
three  longer  cavities.  As  the  small  cavity  turns  in  the  relative  frame  in  the  direction  of  rotation  of 
the  machine  one  qualifies  this  configuration  as  super  synchronous. 

The  consequences  of  these  instabilities  are  dual:  on  the  one  hand,  they  induce  pressure  fluctuations 
upstream  and  downstream  from  the  inducer  and  on  the  other  hand,  they  generate  radial  loads  that  can  be 
responsible  of  the  crushing  of  the  bearings. 

1.2  Two-Phase  Flow  Fluid  Mechanics 

This  paragraph  is  devoted  to  the  presentation  of  the  basics  of  two  phase  flow  fluid  dynamics. 

1.2.1  One  Fluid  Conservation  Equations 

Fluid  dynamics  is  governed  by  mass,  momentum  and  energy  conservation  principles;  these  principle  are 
classically  expressed  by  the  following  local  equations: 

f+v-^)=° 

+  V.(/Pu  ®  U)  =  -S/p  +  V.x  +  pfv 

-  +  V.(phu)  =  -V.q  +  ^-  +  u.S/p  +  </> 
ot  ot 

where  p,  U,  p,  Fv  and  H  are  respectively  density,  velocity,  pressure,  forces  per  unit  of  volume  and 
enthalpy.  These  equations  describe  the  temporal  evolution  of  a  single  fluid  flow.  When  two  fluids  are  put 
together,  an  interface  delimit  the  existence  domain  of  each  phase.  We  examine  next  how  to  model  these 
patterns. 
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1.2.2  Discontinuities 


Let  us  examine  the  effect  of  the  introduction  of  discontinuities  into  the  previous  equations.  We  use  as  a 
convention  [x]  =  x2  -  xt  to  describe  the  variation  of  a  quantity  x  when  crossing  an  interface.  The 
conservation  of  an  extensive  scalar  quantity  per  mass  unit  g  is  classically  written  as: 


D_ 

Dt 


J  pgdv  =  ~\  <P  ■  nclS  +  |  pgvo,dv 

DSD 


pgVr  +  (p  •  n 


—  0  where  Vr  is  the  relative  speed 


This  relation  imposes  when  crossing  an  interface 
to  the  interface. 

•  Mass  conservation: 

The  mass  conservation  corresponds  with  the  previous  notations  to  g=0  and$?=0,  which  is 
written  pA  v \  —vt  •  n  —  p2\  v2  -  vt  •  n  =  0,  that  is  to  say:  mx  —  m2  =  0  where 


mk  -  Pk\^vk  ~  vi  y  n  is  the  mass  flow  of  the  phase  k  through  the  interface.  This  relation 
illustrates  the  conservation  of  the  total  mass  flow  through  the  interface. 

Momentum  conservation: 

The  quantity  to  be  considered  is  g  =  V  and  (p  =  —a  where  cr  is  the  constraints  tensor.  Then, 


one  can  write 
tension  effects. 


pVr  -nV-cr-n 


=  y/cn  +  V  y  .  The  right-side  term  corresponds  to  the  surface 


Energy  conservation: 


The  quantity  to  be  considered  is  g  =  e  +  ^V2  and.  q>  =  q  -  a  •  V  .  Then,  when  neglecting  the 


surface  tension 


effects,  one  obtains:  m[[e  ]  +  1  [ Vr 2  ])  =  - 


q  •  n 


+ 


T  Vr 


Thus,  it  appears  that  the  equations  used  for  traditional  CFD  for  single  phase  flows  can  be  easily  extended 
to  the  diphasic  problems.  However,  if  this  theoretical  extension  remains  possible,  it  also  appears  that 
interfaces  have  to  be  characterised,  which  is  a  harsh  numerical  problem. 

1.2.2  Averages 

Just  like  for  turbulence  modelling,  the  description  of  all  the  diphasic  structures  leads  quickly  to 
multiplicity  of  the  scales  from  temporal  as  well  as  from  spatial  point  of  view.  If  one  is  interested  in  the 
evolution  of  the  macro-scales,  it  is  natural  to  wonder  about  the  relevance  of  an  average  approach.  In 
theory,  three  approaches  can  be  developed:  the  statistical  average,  the  temporal  average  or  the  space 
average.  Just  like  for  turbulence,  the  statistical  average  is  difficult  to  implement  hence  the  recurring  use  of 
the  assumption  of  ergodicity  justifying  the  use  of  temporal  or  spatial  average. 

The  average  approaches  require  the  introduction  of  phase  indicators  Xk  that  are  defined  in  the  following 
way: 
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Xk{x,t) 


fl  if  x  belongs  to  phase  k 
I  0  elsewhere 


This  phase  indicator  function  allows  to  define  the  phase  k  existence  domain  ©  ^  like  the  set  of  points  x 
for  which  Xk(x,t)=  1  at  instant  t.  This  ©  ^  quantity  is  used  to  define  the  volume  fraction  of  phase  k  inside 

^  <Da 

the  whole  domain:  ©  =  /  (Dk  and  ^  .  To  describe  cavitating  flows,  the  more  usual  parameter 

is  void  fraction  a  defined  as  the  vapour  volume  fraction. 

Two  space  averages  are  commonly  used: 

•  The  average  on  the  whole  domain:  (  / )  =  ^  |0  fd<D 

•  The  average  on  the  phase  k  existence  domain:  If)  —  - —  f  fd(D 


From  the  temporal  point  of  view,  a  space  point  x  is  fixed  and  one  defines  the  concept  of  presence  time  of 
the  phase  k  Tk  as  the  set  of  instants  t  when  the  phase  k  is  located  in  x  . 

Two  time  averages  are  commonly  used: 


-  Ip 

The  average  over  the  time  interval  T:  f  —  —  fdt 

T 


•  The  average  over  the  presence  time  of  the  phase  k:  / 


Tk 


Lastly,  by  analogy  with  the  volume  fraction  associated  to  the  spatial  point  of  view,  a  temporal  fraction  can 


be  defined  as  a \ 


Assuming  ergodicity  leads  to  identify  the  averages  and  the  concepts  of  volume 


and  temporal  fraction. 


To  develop  a  numerical  model  with  finite  volume  methods  based  on  an  averaged  approach,  it  is  necessary 
to  apply  the  space  average  operator  to  the  instantaneous  local  equations.  This  operation  is  quite  detailed  in 
many  books  on  two  phase  flows  ([6]  for  example).  These  calculations  lead  for  each  phase  to  the  following 
averaged  equations: 

•  Mass  conservation: 


d^kPk 

8t 


+  ^yXkPkUk 


=  r* 


where  TK  is  the  mass  creation  rate  of  the  phase  k. 

•  Momentum  conservation  for  each  phase: 

8<R-‘^U‘  +  v(r. a ITkFk  )=pkF  +  VTk 
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•  Energy  conservation  for  each  phase: 

dCRj^fk  +  ^pkEkUk^  =  pkF-ul  +  VTk  -UI-  Vqk 


where  Ek  =  Ul  +  ek ,  with  E  energy  internal. 

In  each  equation  some  source  terms  appear  (mass  transfer  Tk  ,  volume  and  surface  forces  FK  TK,  heat  flux 
qK)  that  have  be  modelled:  they  include  in  particular  the  interfacial  effects  averaged  over  the  control 
volume. 


2.0  PHYSICAL  MODELS 

The  models  that  are  developed  to  simulate  cavitation  phenomena  can  be  gathered  in  two  main  classes:  the 
direct  models  and  the  averaged  models. 

2.1  Direct  Diphasic  Models 

We  define  as  direct  diphasic  models,  the  models  that  intend  to  solve  all  the  diphasic  structures.  Two 
approaches  can  be  followed  to  establish  the  equations  to  be  solved:  Lagrange  approach  and  Euler 
approach. 


2.1.1  Lagrange  Approach 

The  Lagrange  approach  consists  in  solving  the  movement  equations  of  a  particle  on  its  trajectory  in  order 
to  calculate  the  evolution  of  its  position  and  speed.  Considering  a  set  of  rigid  particles  for  which  heat  and 
mass  transfer  can  be  neglected,  the  system  to  be  solved  is  the  following: 


d^p 

dt 


=  up 


mp 


du  p  j~i 

~dr=LF‘ 


i, 


deQp  _ rrt 

~dT 


where  mp  is  the  mass  of  the  particle,  Ip  its  moment  of  inertia,  Ft  the  forces  acting  on  the  particle  and  T 
the  torque  generated  by  a  viscous  fluid  on  a  rotating  particle. 

Analytical  solutions  for  the  various  forces  and  the  torque  are  available  only  for  small  particule  Reynolds 
numbers  (Stokes  configuration).  The  extension  of  these  results  to  higher  Reynolds  numbers  is  in  general 
based  on  empirical  correlations. 

The  various  models  can  be  gathered  in  three  classes  according  to  the  level  of  fluid/particle  interactions 
[19].  The  first  level  of  simplification  (one-way  coupling)  consists  in  considering  that  forces  acting  on  a 
particle  are  only  lift  and  drag;  moreover,  the  effect  of  lift  and  drag  is  limited  to  a  single  particule  (no 
action  on  fluid  flow  nor  other  particules).  The  second  level  of  simplification  (two-way  coupling)  consists 
in  adding  the  interaction  with  the  liquid  taking  into  account  mass,  momentum  and  energy  transfers.  Lastly, 
the  most  complex  models  (four-way  coupling)  take  into  account  the  interactions  between  the  particles.  At 
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present  time,  it  seems  that  the  applications  of  this  type  of  methods  to  cavitation  modelling  is  mainly 
limited  to  the  study  of  the  cavitation  nuclei.  The  justification  of  a  level  of  approximation  is  done  according 
to  the  dispersion  of  the  two  phase  flow. 


Distance  inter-particules 


100 


10 


One-way  coupling 


Two-way  coupling 


Four-way  coupling 


4.10" 


4.10" 


Taux  de  vide  a 


Figure  2.1 :  Evolution  of  the  Distance  Inter-Particles  According  to  the  Void  Fraction. 


2.1.2  Euler  Approach 

The  Euler  approach  is  based  on  the  resolution  of  the  Navier-Stokes  equations  combined  with  a 
reconstruction  technique  of  the  interfaces.  For  several  decades,  a  major  effort  has  been  carried  out  on  the 
basis  of  VOF  (Volume  Of  Fluid)  method.  For  the  presentation  of  these  methods,  we  will  place  ourselves 
in  the  case  of  a  two-dimensional  grid  with  of  constant  space  steps  (H  =  Ax  =  Ay).  The  “Volume  Of  Fluid” 

calling  comes  from  the  introduction  of  the  fluid  fraction  C  y  defined  by  C-  h2  =  JJ  H(x ,  y)dxdy  where 

ij 

H  is  the  characteristic  function  of  the  fluid  (H  is  equal  to  1  in  the  fluid  and  0  elsewhere).  Consequently,  C 
takes  values  ranging  between  0  and  1  when  the  mesh  is  cut  by  an  interface,  while  it  equals  0  or  1 
elsewhere.  The  description  of  the  vapour  structures  is  then  relying  on  the  reconstruction  and  propagation 
of  this  interface.  The  two  most  famous  techniques  used  to  determine  the  interface  position  are  SFIC 
(Simple  Line  Interface  Calculation)  and  PLIC  (Piecewise  Linear  Interface  Construction)  methods. 


1 

Figure  2.2:  Interface  Computed  by  SLIC  and  PLIC  Method. 

The  difference  between  the  two  techniques  lies  in  the  fact  that  local  interface  in  each  cell  is  imposed  or  not 
to  be  aligned  with  a  cell  frontier. 

2.2  Averaged  Models 

When  the  void  fraction  values  become  higher  than  4.10"4,  the  interactions  between  particles  must  be  taken 
into  account,  which  generates  a  lot  of  difficulties  to  use  direct  models.  As  the  difficulties  that  have  to  be 
faced  are  similar  to  those  induced  by  turbulence  (multiplication  of  the  space  scales),  we  adopt  the  same 
kind  of  approach  that  has  been  developed  for  turbulence  modelling.  Just  like  for  RANS  turbulence  models, 
we  will  distinguish  several  classes  of  modelling  according  to  the  number  of  equations  used  to  model  the 
diphasic  structures.  The  common  point  of  all  these  models  is  the  use  of  the  mixture  conservation 
equations.  The  variables  in  the  mixture  are  defined  by  averaging  the  variables  of  each  phase: 

•  density:  pm  =  apv  +  (l  -  a)p, 

•  momentum:  ( pu  )m  =  a(pu\  +  (l  -  a\pu\ 
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•  energy:  (pE)m  =  a(pE\  +  (l  -  a\pE\ 

2.2.1  Models  with  n  (n  >2)  Equations 

Having  n  equations  with  n  >2  means  that  the  slip  between  the  phases  (different  speeds  for  the  liquid  and 
the  vapour)  has  to  be  taken  into  account.  The  additional  equations  with  respect  to  the  mixture  conservation 
equations  are  conservation  equations  of  liquid  or  vapour.  The  more  equations  you  try  to  solve  the  more 
closing  terms  you  have  to  model. 


2.2.2  1  Equation  Models 

In  these  models,  a  no  slip  condition  between  the  phases  is  assumed.  The  only  additional  equation  that  has 
to  be  solved  is  the  mass  conservation  of  liquid  or  vapour: 


PM  )  —  0  mass  conservation  of  the  mixture 
+  v(apvu)  =  rv  mass  conservation  of  the  vapour 


The  modelling  effort  then  consists  in  determining  a  formulation  of  the  source  term  rv. 


2.2.2. 7  Use  of  the  Rayleigh-Plesset  Equation 

The  Rayleigh-Plesset  equation  describes  the  evolution  of  the  radius  of  a  spherical  bubble  located  in  an 
infinite  incompressible  fluid  that  is  supposed  to  be  at  steady  state. 

P^UiL  =  R  ™  +  2  (iffi)1  +  4  Yl. M  +  lL_ 

Pl  dt 2  2'dt'  R  dt  Pl  R 


For  practical  modelling  of  cavitation  in  CFD  codes,  this  equation  is  simplified  by  neglecting  surface 
tension,  viscosity  and  second  order  time  dependence: 


dR 

dt 


sgn (Pv  -  P 


Assuming  that  the  bubbles  density  n  is  uniform  and  constant,  the  evolution  of  the  void  fraction  is  straight 
forward  as  the  vapour  volume  is  n^rcR?  .  Concerning  this  approach,  it  has  to  be  kept  in  mind  that  the 


Rayleigh-Plesset  equation  is  not  solved  and  that  the  reference  to  bubbles  is  artificial:  the  most  basic 
expression  of  Rayleigh-Plesset  equation  is  simply  used  to  formulate  a  source  term  for  vaporisation  in  the 
frame  of  a  single  fluid  model. 


2. 2. 2. 2  Empirical  Modelling  of  the  Source  Terms 

To  take  into  account  the  convection  of  void  fraction,  one  can  work  directly  on  modelling  vaporisation  and 
condensation  rate  so  that  they  can  be  included  in  vapour  mass  conservation  equation: 

—ff  +  fptphU  )  =  -[m-  +  m+). 

Following  formulations  proposed  by  Kunz  allowed  to  get  good  numerical  results: 
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•  Vaporisation  rate:  m 


•  Condensation  rate:  m+ 


^jestPj  a min(0,  P~  Pv) 

1  tt2 

-  Pi  Ujx 

CprodP^jl  -  a) 

too 


C  dest  and  C  prod  have  to  be  tuned  by  comparison  to  experimental  data. 


2.2.3  0  Equation  Models 

For  0  equation  models,  only  the  mixture  equations  are  solved.  Consequently,  a  barotropic  state  law 
describing  the  behaviour  of  the  mixture  has  to  be  proposed.  The  density  evolution  is  generally  described 
by  a  piecewise  defined  function: 

r  pL  if  p>p  v  +  Ap 

p=  J  pv  if  p<p  v  -  Ap 

L  P(P)  if  P  v  -  Ap<p<p  v  +  Ap 

where  Ap  is  half-width  of  the  transition  from  vapour  to  liquid.  Three  functions  are  usually  used  to  ensure 
the  transition  between  the  vapour  and  the  liquid:  integration  of  analytical  formula  of  speed  of  sound, 
polynomial  connection  and  sinusoidal  connection. 

•  In  a  liquid  /  incondensable  gas  mixture,  the  speed  of  sound  can  be  easily  calculated  from  the  speed 
of  sound  in  each  phase  and  the  void  fraction  value.  Indeed,  for  a  stationary  diphasic  mixture 
without  mass  transfer: 


a 


-  +  - 


1  -a 


p£  Pici 


,  and  then 


PmCm 

\  =  [apg  +  (l  -  a)p,  ] 


a  +  1  -  a 
PgCg  Pi  cf 


dP 


Assuming  that  the  speed  of  sound  doesn’t  depend  on  entropy,  it  comes  Cm  =  .  Lastly,  the 

dp 

integration  of  this  equation  between  p  v+  Ap  and  the  local  pressure  provides  the  density  value  [7]. 

•  Connection  from  liquid  to  vapour  can  be  defined  with  a  polynomial  connection  [17]: 


j 

p  =  Yj  a-  p‘ 


i= 0 


Coefficients  Ai  are  chosen  to  ensure  continuity  with  liquid  and  vapour  phases;  moreover,  they  are 
selected  in  order  to  obtain  a  sharp  density  gradient  when  cavitation  appears. 

Last,  connection  from  liquid  to  vapour  can  be  defined  with  a  sine  function: 


Pi  +  Pv  Pi  ~  Pv  • 

p  =  — — + 

2  2 


P-Pv 


V  ^*min 


Pi  ~  P 


v  J 
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where  Cmin  corresponds  to  the  minimal  speed  of  sound  in  the  mixture.  This  model  has  been 
developed  at  LEGI  for  CNES  purposes  since  the  early  90’ s  [5].  All  the  results  presented  in  the 
next  paragraphs  are  obtained  with  this  model. 

Two  independent  parameters  are  characterising  this  law:  pv  and  c  min.  The  vapour  density  should  not  be  a 
parameter  since  its  value  is  known.  However,  this  value  is  so  low(0,0256kg/m  3  for  water)  that  it  can 
generate  many  numerical  problems:  the  strong  non  linearity  of  the  law  associated  with  a  very  low  value  pv 
value  requires  a  particularly  strict  numerical  processing.  Consequently,  pv  becomes  a  parameter  of  the  law 
whose  influence  must  be  estimated.  Moreover,  as  pv  is  not  imposed  at  the  vapour  density  value,  it  has  to 
be  considered  as  the  minimum  density  value  of  the  mixture.  When  comparing  to  experimental  results,  this 
hypothesis  is  well  validated  as  maximum  void  fraction  value  is  around  90  %;  indeed,  some  water  drops 
remain  present  inside  the  cavitation  sheets.  Concerning  the  Cmin  value,  the  problem  is  definitely  more 
complex  since  measurement  of  speed  of  sound  in  cavitating  flows  is  not  so  easy.  Consequently,  this  value 
has  been  tuned  by  comparison  to  experimental  data;  in  our  calculations,  the  value  of  0^  is  2,25m/s. 


3.0  PRESENTATION  OF  THE  SOLVER 

In  this  part,  we  give  a  description  of  Euranus  that  is  the  Fine/Turbo™  solver  in  which  the  cavitation 
models  developed  at  LEGI  have  been  integrated.  Fine/Turbo™  is  a  complete  CFD  package  (grid 
generator,  solver  and  post-processing)  developed  by  Numeca  International. 


3.1  Introduction 

The  spatial  discretization  is  based  on  a  finite  volumes  approach;  numerical  scheme  can  be  either  centred  or 
decentred.  The  integration  is  based  on  an  explicit  time-marching  method.  During  each  iteration  a  four 
steps  Runge-Kutta  scheme  is  used. 

We  first  discuss  about  the  way  (local  or  global)  the  pseudo  time  step  is  calculated.  We  will  always  speak 
about  time  step  within  the  meaning  of  pseudo  time  with  respect  to  the  time-marching  algorithm.  The  step 
in  real  time  will  be  explicitly  called  physical  time  step.  The  choice  of  a  local  time  step  means  that  the  time 
step  value  will  be  different  in  each  cell  for  a  single  time  step.  Indeed,  it  is  the  CFL  (Courant  Friedrichs 
Lewy)  number  that  is  kept  constant  to  guaranty  the  stability.  Consequently,  as  the  time  step  is  given  by 

Al-=CFL  wjiere  q  js  ceu  volume  and  ps  the  spectral  radius  (the  details  on  the  calculation  of  the  time 
£2  A 

step  will  be  given  next),  it  appears  that  the  time  step  value  will  depend  on  the  cell  as  Q  and  ps  are  local 
variables.  If  one  uses  a  global  time  step,  then  the  minimum  value  of  the  local  time  steps  will  be  used  for 
each  cell,  inducing  a  lower  convergence  speed. 


3.1.1  System  of  Equations 


The  fluid  dynamics  equations  presented  in  previous  paragraph  can  be  easily  written  in  a  conservative 
vector  form. 


d 

dt 


j 


UdD.  + 


|  FdS  =  J  STdQ 

dQ  Q 


where  U  is  the  vector  of  the  conservative  variables  U  = 


pu 

pv 

ypw) 


and  ST  the  source  term. 
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3.1.2  Space  Discretization  and  Calculation  of  Fluxes 

The  space  discretization  is  based  on  storage  of  the  variables  at  the  centre  of  the  cells.  The  calculation  of 
fluxes  is  carried  out  by  summing  fluxes  on  each  face  of  the  cell.  Convective  fluxes  are  calculated  as 

•  n  -  d  , 


■  /  \  1 

(  ^  \ 

F\u)n  Li  =  F 

Ui  +  Ui+ 1 

\  J  1+2 

2 

^  J 

d  j  is  numerical  dissipation  that  is  necessary  to  stabilize  the  centred  schemes.  Its  calculation  utilizes 


differences  of  the  conservative  variables  at  2nd  and  4th  order: 


f 


,[4] 

'i+1 


\ 


A U  3  -A U  , 

i+  ■  i+- 

V  2  2 


-  £\ 


d  ,  =sll\AU  ,  +  e\ 

i+  i- —  i+  ■ 

2  2  2 

where  A  is  the  centred  differential  operator:  A  U  ,  =  Ui+i  -  U ;  . 


[4] 


A 


A U  ,  -  At/  , 


Z+- 

V  2 


i — 

2  J 


Coefficients  s[ 2 1  ands[4]  are  respectively:  s12^  =  3-  ul,2lA  + ,  max^.j,  st,  si+h  si+2) 


£l4\  =  -  maxi 

'  +  2 


r  .  ~  A 

n  1 ,  >(4]  7  -  rm 

v  ’  2  0  i+i  ei+\ J 


and  s]4]  =  — 


1 


r 

£  ,  , 

«+—  l - 

V  2  2J 


Default  values  of  t>p21and  o,1,4  are  1.0  and  0.1.  The  s,  coefficients  are  activators  proportional  to  the 


pressure  gradients:  v  = 


Pm-  2 Pi  +  Pi- 


Pi  + 1  +  2P,  +  Pi  - 1 


Lastly,  the  variable  A,  is  equal  to  the  spectral  radius  (maximum  eigenvalue)  of  the  convective  fluxes 
matrix  multiplied  by  the  surface  and  corrected  with  a  function  introduced  by  Martinelli  and  Jameson 
(1988)  to  increase  dissipation  in  the  directions  where  the  grid  is  stretched. 


3.1.3  Temporal  Integration  and  Dual  Time  Step 

Dual  time  stepping  consists  in  using  for  the  resolution  of  the  unsteady  equations  an  additional  pseudo-time 
derivative: 


III  wdv+illludv  +  llFds  =  lll STdv 

V  V  S  V 

Thus,  each  physical  time  step  is  regarded  as  a  stationary  problem  whose  solution  is  obtained  via  iterations 
in  pseudo-time.  When  convergence  is  reached,  the  pseudo-time  derivative  is  null  and  the  traditional 
unsteady  equation  is  verified.  In  the  particular  case  of  a  steady  calculation,  the  equation  to  be  solved  does 
not  contain  the  temporal  terms: 

lllwdV  +  \l?Ts=H\STdV 

V  S  V 

Thus,  it  is  noted  that  the  discretized  problem  (steady  or  unsteady)  is  always  locally  expressed  as: 
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Q  +  R(u)  =  0 

or  v  7 

By  definition,  R  is  called  residual.  Temporal  integration  (in  fact,  pseudo  temporal)  is  carried  out  for  each 
iteration  by  an  explicit  4  steps  Runge-Kutta  process. 


3.1.4  Preconditioning 

For  low  Mach  numbers  the  time  marching  algorithms  are  very  slow.  Indeed,  convective  eigenvalues  U  and 
acoustic  eigenvalues  U  ±C  become  very  different  when  the  Mach  number  tends  to  0  (low  compressibility). 
One  of  the  solutions  to  this  problem  is  preconditioning:  it  consists  in  modifying  the  temporal  derivative  in 
the  Navier-Stokes  equations  to  obtain  eigenvalues  centred  around  U.  As  this  term  is  not  physical  (it  is 
cancelled  at  convergence),  the  solution  is  not  modified.  On  the  other  hand,  a  very  significant  increase  in 
convergence  speed  can  be  obtained.  In  the  case  of  Euranus,  the  preconditioning  was  mainly  developed  by 
Hakimi  [8]. 

JJJr  1  f dV + 1  !llUdV + F -Ts = JIf STdV 

v  V  L  Ui  y  S  y 


The  resolution  of  a  stationary  problem  is  carried  out  by  removing  the  temporal  derivative.  The  Q  vector 
corresponds  to  a  new  set  of  non  conservative  variables  and  T'1  is  the  preconditioning  matrix 


Q  = 


c  ^ 

Pi 

u 

V 

\WJ 


and  r  1  = 


p2 

(l  +  cc  )w 

J1 

(l  +  a)v 

J1 

(l  +  a)w 

J1 


pg=  p-p  ref  is  a  gauge  pressure  that  limits  the  rounding  errors. 
Due  to  this  modification  of  the  equations,  eigenvalues  become: 

V  -n(l-a)±^  -rc(l 


V-n  and  4- 


0  0  0 

p  0  0 

0  p  0 

0  0  p 


a)  I  +4 n  P2 


Two  parameters  are  used  in  this  formulation  of  preconditioning:  a  that  is  in  general  equal  to  -1  and  (5  that 
depends  on  a  reference  velocity  ( j]-= ■[ 72  ). 


3.2  Presentation  of  the  Cavitation  Module 

In  this  part,  we  explain  how  the  simulation  of  cavitation  is  integrated  in  the  solver. 

3.2.1  Density  Calculation 

On  the  one  hand,  the  use  of  preconditioning  imposes  to  calculate  the  pressure  from  the  mass  conservation 
equation:  pn+1  =  pn  —  RES  *  At  *  f32 
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In  addition,  the  choice  of  a  barotropic  state  law  allows  to  deduce  directly  the  density  value  from  the 


1  M,i  Pi  T  pv  Pi  pv 

pressure  value:  pn+L  =  — b  sin 


f  pn+l  _  pv 


C 

y  '-'min 


Pi  Pv  J 

In  order  to  smooth  the  density  fluctuations,  two  procedures  are  used:  relaxation  and  dissipation. 


3.2. 1. 1  Relaxation 

When  a  value  of  density  is  calculated  by  application  of  the  barotropic  state  law,  a  relaxation  is  applied: 

pn+ 1  =  pn  +  a(pn+1  -  pn ) 


In  order  to  respect  the  state  law,  the  pressure  is  then  calculated  by  inversing  the  state  law: 

Pn 


f  ,n+l  _  Pi  +  A  ^ 


Pn+l  =  pv  + 


Pi 


Pi  -  A 


Moreover,  an  additional  pressure  relaxation  has  been  added  to  control  the  phenomena  of  vaporization 
(respectively  of  condensation)  starting  from  the  pure  liquid  (respectively  from  the  pure  vapour). 


3. 2. 1.2  Dissipation 

The  second  order  dissipation  is  activated  in  the  zones  of  strong  pressure  gradients.  However,  cavitation  is 
a  phenomenon  which  is  characterized  by  strong  density  gradients.  Consequently,  we  advise  to  modify  the 
calculation  of  second-order  dissipation  by  taking  into  account  density  gradients 


=  ±x 


maxi 


(ep>£p) 


where  sp  =  u,1,2  max(,s'/!b  sf,  sf+l,  sf+2 )  and  sP  =  RMUCAV  max(,v/ 


P  oP  eP  eP  ) 

-1?  9  **i  +  b  **i  +  2  / 


Default  value  of  RMUCAV  is  1.0.  The  sf  coefficients  are  defined  by: 


sf 


Pi+i  ~  2pi  +  Pi- 1 

Pi+ 1  +  2 pi  +  Pi_  i 


3.2.2  Modifications  of  Resolution  Process 

The  transition  from  an  incompressible  flow  to  an  highly  compressible  flow  generates  strong  modifications 
in  the  dynamic  behaviour  of  the  numerical  scheme. 

3.2.2. 1  Eigenvalues 

When  one  considers  the  flow  of  a  compressible  fluid,  the  eigenvalues  of  the  system  are  modified;  that  is 
due  to  the  dependence  of  pressure  to  density  (explicit  dependence  on  speed  of  sound).  Thus,  in  the  two 
phase  areas  one  obtains  the  new  eigenvalues: 
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V-n  and 


V  •  n 


1  -  a  + 


P 


2  A 


+ 


/ 


V-n 

V  v 


1  -  a  + 


p 


2  A 


+  4 


-2 

n 


(v-n] 


The  modification  of  the  calculation  of  the  eigenvalues  changes  the  local  time  step  and  the  second-order 
dissipation  which  are  respectively  inversely  proportional  and  proportional  to  the  spectral  radius. 


3. 2.2.2  Speed  of  Sound 


The  speed  of  sound  appears  explicitly  in  the  expression  of  the  spectral  radius  in  the  two  phase  zones;  it  is 
thus  necessary  to  calculate  it.  Assuming  that  the  speed  of  sound  is  independent  from  entropy  it  comes 


c 


and  using  the  state  law  : 


r  . 

'"min 

Jcos 

'  P-Pv  2  N 

2 

V 

L  Cmin  Pi  ~  Pv  ) 

3. 2.2. 3  Turbulence  Model 


The  calculation  of  turbulent  viscosity  in  the  two  phase  zones  can  be  carried  out  in  two  different  ways. 
Either  one  uses  the  standard  formulas,  or  the  correction  proposed  by  Reboud  [18]: 


Pt 


s  A  +(l -af{pi  ~pv) 


where  jjf  is  standard  turbulent  viscosity. 


4.0  UNSTEADY  TWO  DIMENSIONAL  TEST  CASES 

From  now  on,  we  will  focus  on  numerical  results  obtained  with  the  barotropic  model  of  cavitation  that  has 
been  presented.  We  first  examine  unsteady  results  obtained  in  a  venturi. 

4.1  Computation  Methodology 

This  paragraph  is  devoted  to  the  presentation  of  the  computation  methodology  to  perform  numerical 
simulations  of  unsteady  cavitating  flows.  First,  we  focus  on  boundary  conditions  and  on  the  way  the 
cavitation  is  controlled.  Then,  we  evoke  the  concept  of  convergence  when  performing  unsteady  calculations. 


4.1.1  Boundary  Conditions 

For  incompressible  flows,  boundary  conditions  are  usually  an  imposed  velocity  at  inlet  and  an  imposed 
static  pressure  at  outlet.  With  regard  to  the  boundary  conditions  fixed  on  the  solid  walls,  velocity  is  set  to 
0  (adherence  at  the  wall).  Depending  on  the  turbulence  model  that  is  used  (high  Reynolds  or  low 
Reynolds),  different  grids  have  to  be  built.  Thus,  a  high  Reynolds  model  using  a  wall  function  allows  to 
work  with  grids  for  which  the  first  cell  near  a  wall  is  associated  with  y-i-  values  ranging  between  20  and  50 
while  a  low  Reynolds  model  requires  values  ranging  between  1  and  10. 

4.1.2  Cavitation  Numbers 

P  -  Pva 

The  development  of  cavitation  is  controlled  by  the  non  dimensional  number:  a  =  — - - — .  The 

-2pK 

vapour  pressure  is  a  characteristic  of  the  fluid  and  cavitation  appears  when  the  static  pressure  Pref 
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decreases.  For  our  calculations,  we  use  the  fact  that  the  absolute  value  of  static  pressure  does  not  play  any 
role.  Consequently,  only  the  pressure  gradients  are  effective  so  we  can  fix  a  constant  reference  pressure 
while  varying  the  vapour  pressure.  Practically  static  pressure  is  imposed  at  outlet  and  the  pseudo  vapour- 
pressure  is  modified  during  an  imposed  number  of  physical  time  steps. 


Figure  4.1 :  Evolution  of  the  Vapour  Pressure  during  the  Calculation  Transient. 


With  regard  to  the  value  of  Pvap  that  is  imposed  for  the  first  cavitating  iteration,  a  margin  is  taken  from  the 
appearance  of  vapour.  Thus,  if  AP  is  half- width  of  the  barotropic  law,  the  first  value  of  vapour  pressure  is 
equal  to  Pmim  -2  AP,  where  Pmini  is  the  minimum  value  of  static  pressure  in  the  flow.  Thus,  one  avoids  a 
brutal  appearance  of  vapour  that  could  lead  to  numerical  instability.  This  procedure  highlights  the 
importance  of  using  a  very  well  converged  non  cavitating  calculation  as  an  initial  condition  in  order  to 
avoid  pressure  oscillations  that  would  result  instantaneously  in  the  brutal  appearance  of  vapour.  Moreover, 
the  non  cavitating  result  is  also  used  as  a  reference  to  quantify  the  modifications  induced  by  cavitation. 

4.1.3  Convergence  for  Unsteady  Calculations 

There  is  no  clear  rule  that  ensures  a  strict  convergence  of  an  unsteady  calculation,  nevertheless,  several 
elements  must  be  checked  in  order  to  avoid  coarse  errors: 

•  At  mesh  level,  specific  refinements  must  be  carried  out  in  order  to  be  able  to  describe  the  unsteady 
structures  of  the  flow.  Consequently,  a  simple  refinement  of  boundary  layers  is  not  sufficient. 

•  With  regard  to  the  physical  time  step  value,  it  must  be  a  fraction  of  the  period  of  the  studied 
phenomenon.  In  practice,  we  used  time  steps  equal  to  1/100  of  the  period  estimated  for  the 
phenomenon  to  simulate. 

•  Today,  two  major  criteria  are  used:  number  of  pseudo  time  steps  in  each  physical  time  step  and 
residuals  decrease.  We  will  show  next  that  these  two  methods  are  not  sufficient  to  ensure  a  good 
convergence  of  unsteady  calculations. 

4.2  Venturi  Computations 

A  deep  analysis  of  the  structure  of  the  cavitation  sheets  developing  in  venturi  channels  was  performed  by 
Stutz  [20].  This  work  has  provided  a  great  data  base  for  the  validation  of  our  numerical  work. 

4.2.1  Numerical  Study 

In  this  paragraph,  we  present  the  characteristics  of  the  geometry  and  the  major  numerical  parameters. 
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4.2. 1.1  Global  Presentation 

We  first  present  the  main  characteristics  of  the  geometry: 


Figure  4.2:  Geometry  of  the  Venturi. 

We  performed  a  mesh  convergence  study  on  the  three  following  meshes: 


Table  4.1 :  Sizes  of  the  Grids 


Grid 

number  of  nodes 

1 

81*801 

2 

60*150 

3 

60*265 

A  refinement  of  the  mesh  is  performed  in  the  boundary  layers  and  in  the  throat  section. 
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4.2. 1.2  Influence  of  the  Numerical  Parameters 

The  study  of  the  influence  of  numerical  parameters  was  focused  on  the  determination  of  the  necessary 
conditions  allowing  the  simulation  of  vapour  shedding  downstream  from  a  cavitation  sheet.  First  of  all,  we 
present  a  typical  residual  evolution  for  an  unsteady  calculation. 


Residual 


Physical  time  step 


Figure  4.4:  Convergence  Curve  for  Unsteady  Calculation. 

A  discontinuity  is  observed  at  the  beginning  of  each  time  step  when  residuals  are  initialised;  convergence 
will  be  reached  if  the  number  of  internal  iterations  (pseudo  time  steps)  is  sufficient  to  get  stabilisation  of 
the  residual  at  his  minimum  value. 

Influence  of  the  following  parameters  have  been  analysed: 

•  Transient  duration: 

For  a  number  of  time  steps  higher  than  500  in  the  a  transient,  the  calculation  led  to  steady 
cavitation  at  convergence.  We  were  able  to  exhibit  an  unsteady  behaviour  of  the  cavitation  sheet 
when  imposing  a  transient  of  100  time  steps. 

•  Dissipation  level: 

Influence  of  2nd  order  dissipation  has  been  deeply  analysed.  It  was  found  that  a  too  high  value  of 
numerical  dissipation  leads  to  stabilisation  of  the  computation. 

•  Convergence: 

Convergence  was  validated  by  checking  the  residual  stabilisation  in  each  time  step  and  by 
comparing  the  results  obtained  for  100  and  200  internal  time  steps.  As  no  main  modifications  were 
observed,  it  could  be  concluded  that  the  value  of  100  was  sufficient. 


4. 2. 1.3  Influence  of  the  Boundary  Conditions 

The  main  difficulty  induced  by  the  change  of  the  boundary  conditions  lies  in  the  comparison  of  the  a 
values.  Indeed,  when  the  static  pressure  is  imposed  at  outlet,  a  adownstream  value  based  on  this  quantity  is 
fixed.  When  the  mass  flow  is  imposed  at  outlet,  this  value  of  ad0wnstream  is  no  longer  fixed.  In  order  to 
evaluate  the  importance  of  this  change,  we  carried  out  a  calculation  by  imposing  total  pressure  at  inlet  and 
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mass  flow  at  outlet.  Then,  the  aupstream  value  based  on  inlet  total  pressure  is  kept  constant  in  the  calculation. 
This  inlet  total  pressure  value  was  imposed  in  order  to  get  a  aupstream  value  equal  to  the  mean  value  of  the 
computation  performed  with  static  pressure  imposed  at  outlet.  The  result  of  the  first  computation  was  a 
large  unsteady  cavitation  sheet  whereas  the  result  of  the  second  computation  was  a  little  steady  cavitation 
sheet.  This  result  is  very  important  because  it  shows  that  for  an  unsteady  cavitating  flow,  the  comparison 
with  experimental  results  is  highly  dependant  on  the  way  you  impose  the  a  value. 

4.2.1  A  Influence  of  the  Turbulence  Models 

We  present  now  the  influence  of  the  turbulence  model  on  cavitating  unsteady  calculations.  Computations 
have  been  performed  with  Baldwin-Lomax,  Spallart-Almaras,  Chien  k-s,  Launder-Sharma  k-s,  Yang-Shih 
k-s  and  wall  functions  k-s  models.  Little  discrepancies  in  the  losses  prediction  led  to  high  discrepancy 
when  using  the  absolute  pressure  values  to  calculate  a  values.  The  maximum  variation  has  been  observed 
between  wall  functions  k-s  and  Chien  k-s  models.  The  uncertainty  in  static  pressure  values  led  to  an 
uncertainty  on  aupstream  value  equal  to  0,16,  which  is  huge  when  comparing  to  observed  dependence  of 
cavity  length  to  aupstream  in  the  experiments 


Figure  4.5:  Cavity  Length  vs.  a  in  a  Venture  [20]. 

This  results  highlights  again  the  difficulties  inherent  to  the  rigorous  experimental  validation  of  cavitating 
computations. 

4.2.2  Analysis  of  the  Flow 

We  present  in  this  section  a  brief  analysis  of  the  cavitating  flow  that  was  computed  in  the  venturi.  We 
focus  on  the  throat  area  where  the  cavitation  sheet  appears.  The  study  is  based  on  the  evolution  of  density 
and  velocity  fields. 
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Figure  4.6:  Cycle  of  Vapour  Shedding. 
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Thus,  the  following  mechanisms  appear  to  be  significant  in  the  vapour  shedding  phenomenon: 

•  A  first  phase  is  associated  to  the  decrease  of  the  cavity  length.  This  phase  is  characterized  by  the 
presence  of  a  back  flow  at  the  closure  area  of  the  cavity.  The  inverse  flow  is  localised  near  the 
lower  wall:  this  phenomenon  was  observed  in  experiments  [1].  It  is  from  now  on  admitted  that  it 
plays  a  driving  role  for  the  generation  of  the  vapour  shedding.  This  first  phase  ends  with  the 
reinforcement  of  the  back  flow  and  the  disappearance  of  the  cavity. 

•  A  second  phase  is  characterized  first  by  the  brutal  reattachment  of  the  flow.  Indeed,  the 
disappearance  of  the  cavity  weakens  the  reverse  pressure  gradient  that  was  creating  the  backflow. 
The  vapour  structure  is  then  convected  in  a  vortex  (t=l,0305  s). 

•  The  last  phase  is  characterized  by  the  appearance  and  the  growth  of  a  new  cavitation  sheet  near 
the  throat.  The  cavity  grows  until  the  development  of  a  backflow  in  its  closure  area. 


4.2.3  Convergence  Analysis 

This  last  part  dedicated  to  the  numerical  simulation  of  cavitating  flows  in  a  venturi  is  focused  on  a  simple 
convergence  check  that  has  shown  that  convergence  is  not  so  easy  to  reach. 

The  main  numerical  difficulties  when  performing  cavitating  computations  with  a  preconditioned  coupled 
solver  are  due  to  the  mass  conservation  equation.  Consequently,  it  is  interesting  to  verify  at  each  time  step 
that  mass  conservation  is  well  respected  with  respect  to  global  considerations.  As  an  example,  we  have 
analysed  the  space  evolution  of  mass  flow  in  the  venturi  channel  for  various  time  steps  of  an  unsteady 
cavitating  computation  that  seemed  to  be  converged  (decrease  and  stabilisation  of  the  residuals  for  each 
time  step). 

In  areas  where  the  fluid  has  been  a  pure  liquid  since  more  than  two  physical  time  steps,  the  mass  flow 
value  should  be  constant  (incompressibility  consequence). 


Evolution  du  debit 
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Figure  4.7:  Space  Evolution  of  Mass  Flow  for  Successive  Time  Steps. 
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Each  curve  is  corresponding  to  a  specific  instant.  The  vapour  is  localised  in  the  zones  where  very  sharp 
variations  are  observed.  Thus,  it  appears  that  the  mass  flow  is  not  constant  even  in  areas  where  the  liquid 
is  definitely  incompressible.  The  origin  of  the  problem  was  clearly  identified  as  a  consequence  of  a  too 
low  number  of  internal  iterations  in  each  time  step. 

4.3  HYDROFOIL 

To  finish  this  chapter  concerning  the  numerical  simulation  of  unsteady  cavitating  flows,  we  present  a  work 
that  was  completed  within  the  framework  of  a  discussion  on  cavitation  modelling  for  the  5th  international 
symposium  of  cavitation  CAV2003:  for  this  occasion,  a  numerical  test-case  was  defined  in  order  to 
compare  the  results  obtained  with  various  models. 


4.3.1  Presentation 

The  test-case  is  an  hydrofoil  that  is  placed  in  a  channel  (confined  flow)  and  that  is  defined  by  the 
following  analytical  relation: 


where  the  constant  a0  to  a4  are  respectively: 
a0=  0.11858 
ai  =  -0.02972 
3.2=  0.00593 
a  3  =-0.07272 
a4  =  -0.02207 


The  cord  is  10  cm  ;  concerning  the  channel,  length  is  1  m  and  width  is  40  cm.  Lastly,  it  has  to  be  noted 
that  the  incidence  angle  is  7°.  With  regard  to  the  boundary  conditions,  a  uniform  axial  speed  of  6  m/s  as 
well  as  a  turbulence  intensity  of  0,1%  are  imposed  at  inlet  while  a  uniform  static  pressure  of  10  bars  is 
imposed  at  outlet. 


The  C-mesh  used  for  the  calculations  is  composed  of  45000  cells  (80*571). 
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Figure  4.9:  C-mesh. 

The  first  cell  size  near  a  solid  wall  is  around  50  pm  to  allow  y+  values  on  the  foil  between  4  and  8. 
Computations  are  performed  with  a  standard  k-e  model  with  wall  functions.  The  non  cavitating 
simulations  have  predicted  a  lift  coefficient  of  0,56  and  a  drag  coefficient  of  0,03. 

4.3.2  Analysis  of  Cavitating  Flows  for  a  =  0,8 

We  first  performed  steady  cavitating  calculations  for  a  =  0,8.  Assuming  a  steady  hypothesis,  the  cavitation 
sheet  obtained  at  convergence  was  stable.  The  development  of  this  cavity  induced  a  decrease  in  lift 
coefficient  from  0,056  to  0,046  and  an  increase  in  drag  coefficient  from  0,03  to  0,05.  As  the  cavity  length 
was  close  to  half  the  chord,  unsteady  calculations  were  next  performed.  These  simulations  exhibited  a 
quasi-periodical  behaviour  of  the  flow  that  is  illustrated  by  the  lift  coefficient  evolution. 


Lift  coefficient  vs  time 


Figure  4.10:  Temporal  Evolution  of  the  Lift  Coefficient. 

We  next  present  the  evolution  of  density  field  during  a  cycle. 
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Figure  4.11 :  Time  Evolution  of  Density  Field. 

A  typical  cycle  can  be  divided  in  two  phases:  initially,  the  cavity  grows  while  the  cloud  resulting  from  the 
previous  shedding  is  convected  downstream.  This  first  phase  is  ended  by  the  implosion  of  the  vapour 
structures,  which  initiates  the  second  phase  characterized  by  a  reduction  in  cavity  length.  Strouhal  number 
of  this  periodic  behaviour  is  0,9. 

Lift  coefficient  variations  can  be  explained  by  comparing  the  evolution  of  load  profiles  and  velocity  field 
near  the  trailing  edge  to  visualise  the  deviation  generated  by  the  foil. 


Q. 

o 


Lift  coefficient  vs  time 


Figure  4.12:  Load  Profiles  for  Various  Instants  in  a  Lift  Cycle. 
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It  appears  that  the  minimum  of  lift  coefficient  is  reached  (red  and  violet  curves)  when  a  load  inversion 
occurs  on  second  half  of  the  chord.  These  instants  correspond  to  t=  1,54s  and  t=l,62s  of  figure  4.11.  It  can 
then  be  concluded  that  the  worst  cavitating  configuration  is  not  obtained  for  the  longest  cavity  but  when 
the  vapour  cloud  reaches  the  trailing  edge. 


Figure  4.13:  Velocity  Field  at  Trailing  Edge. 

Black  lines  are  the  frontiers  of  two  phase  areas.  This  unsteady  cavitating  configuration  led  to  a  mean  lift 
coefficient  of  0,46  and  a  mean  drag  coefficient  of  0,078.  Consequently,  it  appears  that  the  lift  computed  in 
the  steady  calculation  is  a  good  approximation  of  the  mean  value  obtained  in  the  unsteady  one. 
Concerning  the  drag  prediction,  the  conclusion  is  different  as  the  mean  drag  value  is  56%  higher  in  the 
unsteady  computation  than  the  value  obtained  in  the  steady  calculation. 


5.0  NUMERICAL  METHODOLOGY  FOR  TURBOMACHINERY 

In  this  chapter,  we  focus  on  numerical  specificities  of  the  cavitating  simulations  in  turbomachinery. 

5.1  Presentation 

Two  main  topics  are  discussed  in  this  first  part.  First,  the  limits  of  the  classical  convergence  criteria  are 
discussed.  Then,  the  influence  of  the  computational  methodology  is  examined. 

5.1.1  Discussion  on  Convergence 

From  a  theoretical  point  of  view,  convergence  is  well  defined  in  the  following  way:  “the  discrete 
approximation  is  convergent  if  the  unique  discrete  solution  tends  towards  the  exact  solution  when 
decreasing  the  time  step  At  and  the  space  step  Ax  towards  0.”  In  practice,  it  is  impossible  to  verify  these 
conditions  when  numerical  schemes  become  complex  (nominal  case  for  CFD  computations). 
Consequently,  to  characterise  the  convergence  level  of  a  calculation,  some  practical  criteria  have  been 
defined  that  are  based  on  the  evolution  of  the  residuals.  Indeed,  these  variables  represent  exactly  the  error 
that  is  made  in  the  conservation  equations.  The  more  common  criterion  is  related  to  the  decrease  of  the 
residuals.  Usually  if  the  residuals  decrease  of  4  orders  of  magnitude  then  calculation  may  be  considered  as 
converged.  This  evaluation  is  carried  out  by  defining  a  global  residual  from  the  local  residuals: 
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The  operator  moyRMs  corresponds  to  RMS  average  and  each  residual  RES  is  normalised  with  its  initial 

value.  Consequently,  one  will  consider  that  the  residuals  have  decreased  of  4  decades  when  RESrms  =  -4. 
We  think  that  this  criterion  requires  a  careful  application  for  two  reasons: 

•  On  the  one  hand,  the  criterion  is  depending  on  the  initial  conditions:  a  badly  initialized  calculation 
will  have  a  great  decrease  of  the  residuals  without  ensuring  a  good  convergence. 

•  On  the  other  hand,  the  criterion  does  not  take  into  account  the  dynamic  evolution  of  the  residuals. 
Thus,  this  criterion  does  not  consider  as  converged  a  calculation  even  if  the  solution  doesn’t 
evolve  any  more. 

5.1.2  Methodology  Impact  on  Head  Drop  Results 

We  present  in  this  paragraph  two  major  aspects  of  computational  methodology  for  the  simulation  of 
cavitating  flows: 

•  The  choice  of  a  procedure:  step  by  step  or  continuous. 

•  The  choice  of  a  formulation  of  the  equations:  steady  or  unsteady. 

In  the  past,  head  drop  simulations  were  performed  by  decreasing  continuously  the  cavitation  number  and 
by  using  unsteady  formulation  of  the  equations.  We  advise  not  to  use  this  kind  of  computations  but  to 
perform  multiple  steady  computations.  The  main  objective  is  to  master  the  convergence  level  by  using 
robust  criteria.  To  confirm  this  statement,  we  examine  some  of  the  consequences  of  an  unsteady 
continuous  procedure: 

•  The  main  reason  explaining  the  use  of  unsteady  formulation  is  a  lack  of  stability.  Indeed,  the 
temporal  derivative  acts  as  a  compensation  term  with  respect  to  the  brutal  variations  of  pressure 
and  density. 

•  A  continuous  variation  of  vapour  pressure  generates  an  anticipated  prediction  of  head  drop.  As  a 
proof  element,  we  present  the  typical  evolution  of  pressure  rise  for  a  steady  cavitating 
computation  of  a  rocket  engine  inducer  (inducer  2  presented  in  §7). 


Figure  5.1 :  Evolution  of  Total  Pressure  Rise  in  the  Course  of  Calculation. 


3-26 


RTO-EN-AVT-1 43 


£ 

NATO 

31 

OTAN 

Numerical  Modelling  of  Cavitation 


Thus,  it  appears  clearly  that  pressure  rise  decreases  first  when  cavitation  is  developing  (the  transient  on  the 
vapour  pressure  is  performed  during  the  first  5000  iterations).  When  the  decrease  in  vapour  pressure  is 
stopped,  a  pressure  rise  increase  is  observed,  which  proves  that  convergence  was  far  from  being  reached  at 
each  time  step. 

5.2  Problems 

Simulation  of  cavitating  flows  in  turbomachinery  generates  specific  problems  with  respect  to  numerical 
stability.  The  major  problem  is  the  transition  from  pure  liquid  to  pure  vapour  of  a  cell  in  an  single  time 
step  that  can  provoke  brutal  divergences.  We  next  present  how  that  occurs: 

•  We  first  mention  that  density  is  calculated  from  the  pressure  ;  then,  an  abrupt  variation  of  density 
is  caused  by  a  pressure  step.  As  the  pressure  is  determined  by  mass  conservation  equation,  one 
deduces  that  the  dynamic  behaviour  of  this  equation  is  responsible  for  the  divergence. 

•  The  mass  conservation  equation  manages  the  evolution  of  the  pressure  by  the  following  process: 
pn+l=pn-RES*AT*/3 2 .  RES  corresponds  to  mass  flow  balance  in  the  cell  and  Ax  is  the  pseudo 

time  step.  Then,  this  equation  converts  instantaneously  an  error  in  mass  conservation  into  a 
pressure  variation. 

•  Lastly,  instability  is  due  to  sharp  variations  of  density,  which  means  that  the  problem  can  be 
identified  as  local  with  respect  to  the  pressure  field.  Indeed,  because  of  the  stiffness  of  the  state 
law,  the  transition  from  pure  liquid  to  pure  vapour  can  be  caused  by  a  variation  of  static  pressure 
of  about  0,1  bar  what  is  low  compared  to  the  pressure  variations  observed  between  two  iterations. 

5.3  Solutions 

In  this  paragraph,  we  detail  the  various  ways  that  were  studied  to  improve  the  code  stability.  The  P 
coefficient  is  a  fixed  parameter  calculated  using  the  following  formula: 

P2=Pfflref  where  /?02=  3  and  Uref  is  a  reference  velocity 

Consequently,  we  initially  tried  to  decrease  the  reference  velocity  value.  Generally,  it  is  recommended  to 
use  the  highest  velocity  of  the  flow  in  order  to  maximize  artificial  dissipation.  Initially,  we  decided  not  to 
modify  the  calculation  of  the  spectral  radius  in  order  to  preserve  high  levels  of  dissipation.  Stability  was 
indeed  improved  but  the  results  obtained  were  characterized  by  very  bad  convergence  level  (2%  of  error 
on  the  mass  flow  conservation)  and  by  cavity  shapes  far  from  being  probable. 


Figure  5.2:  influence  Dissipation  on  the  Cavity  Shape. 
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The  increase  in  dissipation  generates  a  reduction  of  cavity  length  associated  with  a  very  strong  smoothing 
of  the  interface.  Consequently,  we  decided  to  modify  the  calculation  of  the  spectral  radius  in  a  coherent 
way  with  the  reduction  of  (3  and  a  great  improvement  was  observed.  Then,  in  order  to  minimize  the 
computing  time,  we  decided  to  limit  the  influence  of  the  P  reduction  in  cavitating  areas.  Finally,  we 
defined  a  P  law  depending  on  static  pressure. 


Figure  5.3:  Evolution  of  p  versus  Static  Pressure. 


Linking  between  the  minimum  and  maximum  value  of  P  is  made  outside  from  the  barotropic  state  law  in 
order  to  stabilize  the  evolution  of  the  cells  where  pure  liquid  is  converted  into  vapour.  This  law  has 
provided  high  stability  to  the  algorithm. 

5.4  Influence  of  State  Law  Parameters 

We  are  mainly  interested  in  the  influence  of  the  parameters  of  the  state  law  :  RHOV  and  AMIN.  It  is 
important  to  keep  in  mind  that  RHOV  should  not  be  considered  as  the  vapour  density  value  but  rather  as 
the  minimal  density  value  of  the  mixture;  thus,  a  value  of  100  for  cold  water  means  that  the  maximum 
void  fraction  value  that  we  can  simulate  is  90  %,  which  is  higher  than  the  maximum  void  fraction  values 
measured  in  a  venturi. 
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The  value  of  500  allows  to  predict  a  head  drop  whose  uncertainty  is  about  10  %.  On  the  other  hand,  the 
operating  point  calculated  with  RHOV=50  is  very  close  to  that  calculated  for  RHOV=100  (deviation  lower 
than  2%). 

The  study  of  the  AMIN  influence  was  also  led  on  this  geometry  for  values  going  from  30  to  100.  No  effect 
could  be  observed  on  the  total  performance  of  the  impeller.  However,  we  observed  a  significant 
dependence  of  the  cavity  volume  on  this  parameter.  For  example,  a  variation  of  AMIN  from  35  to  75  led 
to  an  increase  of  cavity  thickness  of  about  20  %. 


6.0  HEAD  DROP  OF  IMPELLERS 

SHF  impeller  was  designed  by  NEYRPIC  company  in  the  frame  of  a  research  project  on  partial  flow  rates 
operation.  This  impeller  was  tested  at  ENSAM  Lille  in  air,  at  INSA  Lyon,  EPF  Lausanne  and 
HYDRO  ART  Milan  in  water.  Numerical  studies  were  also  performed  by  Combes  [2].  The  main 
characteristics  of  SHF  impeller  are: 

•  A  number  of  blades:  7. 

•  Number  of  revolutions:  3000  tours/min. 

•  Nominal  capacity:  0,157  m3/s. 

•  Diameter  piping  upstream:  200mm. 

•  Diameter  left:  330  mm. 
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6.1  Mesh  and  Cutting  Plans 

In  this  paragraph,  we  quickly  recall  the  main  characteristics  of  the  mesh  that  is  used  for  the  calculations. 

•  Size  of  first  cell  near  a  wall  =10  pm. 

•  Additional  refinement  at  leading  edge  with  an  expansion  coefficient  of  0.3. 

•  10  cells  without  size  variations  in  the  center  of  the  channel  in  azimuth  direction. 

•  12  cells  without  size  variations  in  the  center  of  the  channel  in  radial  direction. 

Total  number  of  internal  nodes  is  250000  ;  this  value  has  been  chosen  after  a  mesh  convergence  based  on 
computations  with  other  meshes  of  170  000,  500000  and  1000000  nodes.  The  boundary  conditions  are: 
imposed  velocity  at  inlet,  imposed  static  pressure  at  outlet,  and  periodicity  conditions. 


Figure  6.2:  3d  and  Meridian  View  of  the  Mesh. 


To  study  the  local  flow  structures,  we  created  cutting  plans  orthogonal  to  the  meridian  direction.  All  the 
results  will  be  presented  with  the  following  conventions: 
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Figure  6.3:  Cutting  Plans. 

The  plans  are  numbered  from  I  to  VIII  from  inlet  to  outlet.  Shape  of  the  plans  is  very  distorted  due  to  high 
twisting  of  the  blades. 

6.2  Cavitating  Results 

The  presentation  of  the  cavitating  results  is  carried  out  in  two  steps:  first  of  all,  we  explain  the  head  drop 
phenomenon  from  a  global  point  of  view  and  then,  we  expose  some  local  modifications  induced  by  the 
development  of  the  cavitation  sheets. 


6.2.1  Global  Results 

First,  we  expose  the  elements  explaining  the  head  drop  phenomenon.  For  that,  we  first  present  the  head 
drop  curves  that  we  obtained  for  various  flow  rates: 


Head  drop  at  various  flow  rates 


Figure  6.4:  Head  Drops  for  Various  Mass  Flows. 
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Head  is  here  defined  as  the  static  pressure  rise  between  inlet  and  outlet  sections.  The  last  converged  point 
corresponds  to  a  10  %  head  drop.  Moreover,  one  will  note  that  the  calculated  value  of  NPSH 
corresponding  to  head  drop  is  constant  four  different  calculations,  which  is  a  sign  of  good  convergence. 
Last,  the  good  agreement  with  measured  value  has  to  be  noticed.  Next,  we  examine  the  evolution  of 
cavitation  for  decreasing  a  values. 


Figure  6.5:  Evolution  of  the  Cavities. 

The  cavity  appears  initially  on  suction  side  near  the  hub.  It  develops  mainly  in  radial  and  meridian 
directions  until  reaching  the  throat,  which  provokes  the  brutal  performance  drop.  Lastly,  cavitation 
develops  on  pressure  side  very  close  to  the  shroud;  then  an  important  azimuth  and  radial  extension  is 
observed.  The  two  cavities  (suction  face  and  shroud)  are  very  thin. 

With  regard  to  the  head  drop  of  the  machine,  it  is  interesting  to  dissociate  the  energy  transfer  and  the 
losses.  Consequently,  we  present  the  respective  evolutions  of  the  hydraulic  and  mechanical  power. 
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Figure  6.6:  Hydraulic  and  Mechanical  Power  Evolution. 
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Thus,  it  appears  that  hydraulic  and  mechanical  power  decrease  almost  simultaneously,  which  means  that 
the  decrease  in  pressure  rise  is  mainly  due  to  a  reduction  of  the  torque.  Loss  of  efficiency  can  be  seen  as  a 
second  order  as  the  decrease  is  about  2  %  for  the  last  cavitating  point,  which  is  negligible  in  front  of  the 
torque  reduction.  In  order  to  have  a  better  understanding  of  this  evolution,  we  compare  the  blade  load  at 
mid  span  for  various  values  of  NPSH. 


Figure  6.7:  Load  Evolution  at  Mid  Span. 

The  development  of  cavitation  corresponds  to  the  evolution  from  blue  to  red.  Thus,  it  is  noted  that  the 
torque  reduction  is  connected  very  simply  to  the  growth  of  the  cavity.  Static  pressure  is  limited  to  vapour 
pressure  whereas  one  would  need  lower  pressure  to  preserve  the  same  blade  load.  Moreover,  these  curves 
highlight  an  increase  in  the  load  downstream  from  the  cavity. 

Next,  it  is  interesting  to  examine  how  these  modifications  of  the  blade  loading  change  the  energy  transfers 
inside  the  machine. 
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Figure  6.8:  Evolution  of  Internal  Pressure  Rises. 
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Each  histogram  element  represents  the  total  pressure  rise  of  the  impeller  for  a  particular  calculation. 
Calculation  1  is  calculation  without  cavitation  and  calculation  7  is  the  calculation  carried  out  with  the 
lowest  sigma  value.  Total  pressure  rise  is  then  distributed  in  7  portions  corresponding  to  the  7  intervals 
between  the  cutting  plans  defined  on  figure  6.3. 

This  representation  highlights  the  meridian  evolution  of  total  pressure  rise: 

•  In  the  first  portion  (violet  colour),  one  observes  for  the  first  three  operating  points  an  increase  in 
the  total  pressure  rise  that  causes  a  global  performance  increase.  This  phenomenon  can  be 
explained  by  an  incidence  effect:  the  development  of  the  cavity  thickens  the  suction  side  profile 
inducing  a  local  increase  in  the  incidence. 

•  A  compensation  phenomenon  is  observed  downstream  from  the  cavitation  sheet.  This 
phenomenon  is  particularly  visible  in  portion  4  (light  blue)  for  calculations  4  to  7. 

6.2.2  Local  Results 

In  this  part,  we  are  interested  in  the  modifications  of  the  flow  induced  by  the  development  of  a  cavitation 
sheet.  We  use  the  cuts  defined  on  figure  6.3.  Our  study  is  limited  to  the  comparison  of  non  dimensional 
relative  helicity  for  the  most  developed  cavitating  configuration  (NPSH  «  6,3m). 


Figure  6.9:  Non  Dimensional  Relative  Helicity  in  Plans  I  and  II  (Non  Cavitating). 

In  plan  I,  the  flow  seems  very  heterogeneous  and  it  is  difficult  to  give  an  interpretation  of  the  structures 
observed.  In  plan  II,  one  observes  the  development  of  the  first  secondary  flows.  Thus,  the  effect  of 
meridian  curvature  is  particularly  visible  near  the  suction  side.  On  the  pressure  side,  the  swirling  structure 
develops  with  difficulty  and  remains  confined  near  the  shroud. 


Figure  6.10:  Non  Dimensional  Relative  Helicity  in  Plans  I  and  II  (Cavitating). 


3-34 


RTO-EN-AVT-1 43 


1 u 

NATO 

31 

OTAN 

Numerical  Modelling  of  Cavitation 


The  observation  of  the  non  dimensional  relative  helicity  in  plan  I  illustrates  a  very  significant  modification 
of  the  secondary  flows.  In  non  cavitating  configuration,  it  was  difficult  to  distinguish  organized  structures 
whereas  the  presence  of  the  cavity  generates  a  swirling  zone  near  the  interface.  In  plan  II,  one  notes  that 
the  structure  observed  in  plan  I  almost  disappeared  and  that  it  persists  only  near  the  hub.  On  the  other 
hand,  one  notes  the  development  of  a  major  swirling  structure  near  the  suction  side  whose  presence  may 
be  explained  by  the  centrifugation  of  the  wake  of  the  cavity. 


7.0  HEAD  DROP  OF  INDUCERS 

This  part  is  devoted  to  the  simulation  of  cavitating  flows  in  two  rocket  engine  inducers. 


Figure  7.1 :  inducer  Geometries  (left  =  Inducer  1  and  right  =  Inducer  2). 

Both  inducers  were  designed  for  research  objectives  ;  the  first  inducer  was  designed  by  Snecma  and  the 
second  inducer  was  designed  by  Pr.  Tsujimoto  team  from  Osaka  University. 

7.1  Inducer  1 

7.1.1  Mesh  and  Cutting  Plans 

In  this  paragraph,  we  quickly  recall  the  main  characteristics  of  the  mesh  that  is  used  for  the  calculations. 

•  Size  of  first  cell  near  a  wall  =  5  pm. 

•  41  cells  in  azimuth  direction. 

•  41  cells  in  radial  direction. 

•  203  cells  in  axial  direction. 

Total  number  of  internal  nodes  is  350000  for  one  blade  to  blade  channel;  this  value  has  been  chosen  after  a 
mesh  convergence  based  on  computations  with  other  meshes  up  to  1000000  nodes.  The  boundary 
conditions  are:  imposed  velocity  at  the  mesh  inlet,  imposed  static  pressure  at  the  outlet,  and  periodicity 
conditions.  The  turbulence  models  is  standard  k-s  with  wall  functions.  Tip  gap  has  been  taken  into  account 
only  for  mesh  convergence  in  non  cavitating  conditions.  Numerical  results  were  validated  by  comparison 
with  experimental  results  on  static  pressure  evolution  on  the  shroud  for  non  cavitating  conditions. 
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Figure  7.2:  Evolution  of  Static  Pressure  vs  Axial  Length. 


To  study  the  local  flow  structures,  we  created  cutting  plans  orthogonal  to  the  meridian  direction.  All  the 
results  will  be  presented  with  the  following  conventions: 


Figure  7.3:  Cutting  Plans. 


7.1.2  Cavitating  Results 

The  presentation  of  the  cavitating  results  is  carried  out  in  two  steps:  first  of  all,  we  explain  the 
phenomenon  of  head  drop  from  a  global  point  of  view  and  then,  we  expose  the  local  modifications 
induced  by  the  development  of  the  cavitation  sheets. 

7. 7.2. 7  Global  Analysis 

First,  we  expose  elements  related  to  head  drop  phenomenon. 
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Figure  7.4:  Head  Drop  at  Nominal  Flow  Rate. 

Concerning  the  prediction  of  sigma  value  associated  with  head  drop,  a  good  agreement  between  numerical 
and  experimental  values  is  obtained.  The  brutal  behaviour  of  the  head  drop  is  also  well  described.  The 
absolute  head  gap  can  be  explained  by  our  neglecting  the  tip  gap  in  these  calculations. 

Next,  we  examine  the  evolution  of  cavitation  for  decreasing  a  values. 


Figure  7.5:  Evolution  of  the  Cavity. 

Cavitation  appears  initially  near  the  shroud  on  the  suction  side  and  remains  very  localised  for  a  wide  range 
of  a  values.  Then,  cavitation  mainly  develops  in  the  back  flow  creating  a  highly  three-dimensional 
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structure  near  the  shroud.  One  then  attends  the  appearance  of  a  cavity  attached  to  the  leading  edge  along 
the  span  of  the  blade.  This  cavity  grows  mainly  in  the  meridian  direction.  It  will  be  noted  in  addition  that 
the  thickness  of  the  cavities  is  low  and  that  spanwise  extension  remains  rather  uniform.  Then,  we  present 
the  evolution  of  hydraulic  and  mechanical  power. 


Figure  7.6:  Evolution  of  Hydraulic  and  Mechanical  Power. 

It  appears  that  the  drops  of  hydraulic  and  mechanical  power  are  almost  simultaneous,  which  means  that 
the  decrease  in  pressure  rise  is  mainly  due  to  a  reduction  of  the  torque  like  for  the  SHF  impeller.  The  fall 
of  efficiency  is  about  3  %  for  the  last  calculated  point. 

In  order  to  have  a  better  understanding  of  this  evolution,  we  compare  the  blade  load  at  mid  span  for 
various  sigma  values. 


Evolution  de  la  charge 


Corde  normalisee 


Figure  7.7:  Evolution  of  the  Blade  Load. 
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The  development  of  cavitation  corresponds  to  the  evolution  from  blue  to  red.  Thus,  it  is  noted  that  the 
torque  reduction  is  connected  as  for  the  SHF  impeller  to  the  growth  of  cavity  length.  Static  pressure  is 
limited  to  vapour  pressure  whereas  lower  pressures  are  needed  to  preserve  the  same  blade  load.  On  the 
other  hand,  the  effects  of  compensation  observed  on  the  load  profiles  for  SHF  impeller  are  not  obvious 
here.  Moreover,  contrarily  to  what  was  observed  on  SHF  impeller,  one  notes  a  significant  modification  of 
the  static  pressure  profile  on  pressure  side  that  is  responsible  for  part  of  the  head  drop.  Next,  it  is 
interesting  to  examine  how  these  modifications  of  the  blade  loading  change  the  energy  transfers  inside  the 
machine. 
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Figure  7.8:  Evolution  of  Internal  Pressure  Rises. 


Each  element  of  the  histogram  represents  the  total  pressure  rise  of  the  inducer  for  a  particular  calculation. 
Calculation  1  is  performed  in  non  cavitating  conditions  and  calculation  9  is  the  calculation  carried  out  with 
the  lowest  sigma  value.  Total  pressure  rise  is  then  distributed  in  7  portions  corresponding  to  the  7  intervals 
between  the  cutting  plans  defined  on  figure  7.3.  This  representation  highlights  the  meridian  evolution  of 
total  pressure  rise: 

•  In  the  first  portion  (violet  colour),  the  pressure  rise  increase  observed  on  the  SHF  impeller  for  the 
highest  sigma  values  doesn’t  appear. 

•  In  the  second  portion,  one  observes  for  the  second  point  a  decrease  of  pressure  rise  that  is 
responsible  of  the  global  pressure  rise  decrease  observed  on  the  inducer  at  this  point.  Indeed,  it  is 
noted  that  the  pressure  rise  provided  by  other  portions  is  not  affected. 

Finally,  as  it  was  said  by  examining  the  load  profiles,  this  study  confirms  that  the  compensation 
phenomena  observed  on  SHF  impeller  don’t  exist  in  this  inducer. 

7. 7.2.2  Local  Analysis 

In  this  part,  we  focus  on  the  modifications  of  the  flow  induced  by  the  development  of  a  cavitation  sheet. 
We  use  the  cutting  plans  presented  on  figure  7.3  to  visualise  axial  velocity  and  non  dimensional  relative 
helicity  distributions. 
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•  Axial  velocity  : 


Figure  7.9:  Axial  Velocity  in  Plan  I  (left  =  Non  Cavitating,  right  =  Cavitating). 

In  plan  I,  one  notes  that  the  presence  of  the  cavity  generates  heterogeneities  in  the  axial  velocity  distribution. 
In  particular,  one  will  note  the  presence  of  low  velocity  zones  on  suction  side  that  are  due  to  the  obstruction 
generated  by  the  cavity.  In  plan  II,  we  could  observe  that  the  flow  remains  very  homogeneous. 

•  Secondary  flows: 


Figure  7.10:  Distribution  of  Relative  Helicity  in  Plans  I  and  II  (Cavitating). 

In  plan  I,  the  distribution  of  non  dimensional  relative  helicity  is  very  disturbed  on  suction  side  near  the 
shroud.  This  phenomenon  can  be  explained  by  the  presence  in  this  zone  of  the  wake  of  the  cavity.  Near  the 
hub,  the  secondary  flows  are  close  to  those  observed  in  non  cavitating  calculations. 


Figure  7.11 :  Distribution  of  Relative  Helicity  in  Plan  IV  (left  =  Non  Cavitating,  right  =  Cavitating). 
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In  plans  II  and  III,  we  noted  that  the  pressure  side  and  suction  side  vortices  were  very  close  to  those 
observed  in  non  cavitating  computations.  In  plan  IV,  significant  differences  have  been  seen  with  non 
cavitating  distribution.  For  example,  we  could  observe  that  the  pressure  side  vortex  deeply  modified  the 
suction  side  flow  near  the  trailing  edge.  Indeed,  one  notes  disappearance  on  this  face  of  the  negative 
helicity  structure  that  was  convected  towards  the  centre  of  the  channel.  Lastly,  one  also  raises  the 
existence  of  a  small  structure  near  the  hub  that  didn’t  exist  in  non  cavitating  computation.  Consequently,  it 
has  to  be  reminded  that  the  outlet  flow  is  highly  disturbed  by  the  development  of  the  cavity. 

7.2  Inducer  2 

7.2.1  Mesh  and  Cutting  Plans 

In  this  paragraph,  we  quickly  recall  the  main  characteristics  of  the  mesh  that  is  used  for  the  calculations. 

•  Size  of  first  cell  near  a  wall  =10  pm. 

•  41  cells  in  azimuth  direction. 

•  41  cells  in  radial  direction. 

•  149  cells  in  axial  direction. 

Total  number  of  internal  nodes  is  250000  for  one  blade  to  blade  channel;  this  value  has  been  chosen  after  a 
mesh  convergence  based  on  computations  with  other  meshes  whose  cells  number  was  up  to  500000.  The 
boundary  conditions  are:  imposed  velocity  at  inlet,  imposed  static  pressure  at  outlet,  and  periodicity 
conditions.  The  turbulence  model  is  standard  k-s  with  wall  functions.  Tip  gap  has  not  been  taken  into 
account. 


Figure  7.12:  Geometry  and  Definition  of  Cutting  Plans. 


7.2.2  Cavitating  Results 

The  presentation  of  the  cavitating  results  is  focused  on  the  head  drop  phenomenon.  We  first  recall  the  head 
drop  curve  that  was  obtained  at  nominal  flow  rate. 


RTO-EN-AVT-1 43 


3-41 


Numerical  Modelling  of  Cavitation 


Chute  de  performance 


sigma 


Figure  7.13:  Head  Drop. 


We  then  examine  the  evolution  of  cavitation  for  various  sigma  values. 


Figure  7.14:  Evolution  of  the  Cavities. 


Cavitation  appears  on  all  the  height  of  the  blade  with  a  triangular  shape  whose  base  is  on  the  shroud.  This 
cavity  grows  uniformly  when  decreasing  the  a  value.  Cavitation  develops  in  the  back  flow,  which  creates 
a  highly  three-dimensional  diphasic  structure  near  the  shroud  and  extending  upstream.  Cavitation  sheets 
are  thicker  than  for  inducer  1. 
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To  characterize  the  performance  drop,  we  examine  the  evolution  of  hydraulic  and  mechanical  power. 


Evolution  des  puissances  regue  et  fournie 


sigma 


Figure  7.15:  Evolution  of  the  Hydraulic  and  Mechanical  Power. 

It  appears  that  the  hydraulic  and  mechanical  power  drops  are  almost  simultaneous,  which  means  that  the 
decrease  in  pressure  rise  is  mainly  due  to  a  reduction  of  the  torque  like  for  the  SHF  impeller.  The 
efficiency  loss  is  about  3  %  for  the  last  calculated  point.  These  results  are  identical  to  what  has  been  seen 
on  inducer  1 . 

In  order  to  have  a  better  understanding  of  this  evolution,  we  compare  the  blade  load  at  mid  span  for 
various  sigma  values. 


Figure  7.16:  Modification  of  the  Blade  Load  at  Mid  Span. 
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The  development  of  cavitation  corresponds  to  the  evolution  from  blue  to  red.  Thus,  it  is  noted  that  the 
torque  reduction  is  connected  as  for  SHF  impeller  and  inducer  1  to  the  growth  of  cavity  length.  Static 
pressure  is  limited  to  vapour  pressure  whereas  lower  pressures  are  needed  to  preserve  the  same  blade  load. 
The  compensation  effects  observed  downstream  from  the  cavity  for  SHF  impeller  are  lower  on  inducer  2 
than  on  SHF  impeller.  This  effect  didn’t  exist  on  inducer  1;  this  difference  is  explained  by  the  fact  that 
pressure  side  profile  isn’t  modified  contrarily  to  what  was  observed  on  inducer  1.  This  evolution  of  the 
load  also  explains  why  the  head  drop  of  the  inducer  1  is  so  progressive:  this  inducer  was  dimensioned  with 
high  incidence  angles,  which  results  in  high  loading  of  the  first  part  of  the  blade.  Consequently, 
appearance  of  cavitation  in  this  area  generates  a  torque  drop  that  cannot  be  compensated.  Next,  it  is 
interesting  to  examine  how  these  modifications  of  the  blade  loading  change  the  energy  transfers  inside  the 
machine. 


Figure  7.17:  Evolution  of  Internal  Pressure  Rises. 

Each  element  of  the  histogram  represents  the  total  pressure  rise  of  the  inducer  for  a  particular  calculation. 
Calculation  1  is  calculation  without  cavitation  and  calculation  8  is  the  calculation  carried  out  for  the 
lowest  sigma  value.  Total  pressure  rise  is  then  distributed  in  5  portions  corresponding  to  the  5  intervals 
between  the  cutting  plans  defined  on  figure  7.12.  This  representation  highlights  the  axial  evolution  of  total 
pressure  rise: 

•  The  incidence  loading  is  clearly  highlighted  on  the  energy  transfers  histograms  since  one  can  note 
that  the  transfers  taking  place  in  the  last  three  sections  are  weak  in  front  of  those  in  the  first  two 
sections. 

•  In  the  first  section,  a  pressure  rise  increase  is  observed  while  cavitation  is  developing.  This  is 
similar  to  what  was  found  on  SHF  impeller. 

•  The  compensation  phenomenon  observed  downstream  from  the  cavity  on  the  load  profiles  is  very 
clear  when  comparing  the  evolution  of  the  first  two  sections. 

8.0  INSTABILITIES  IN  A  CASCADE 

This  last  part  of  the  note  is  devoted  to  the  presentation  of  some  numerical  simulation  of  cavitation 
instabilities  in  a  2D  blade  cascade. 
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8.1  Presentation 

To  obtain  the  2D  geometry,  a  3D  inducer  is  cut  at  a  constant  radius  Rc  equal  to  70%  of  the  tip  radius  R. 
The  boundary  conditions  are:  imposed  velocity  at  inlet,  imposed  static  pressure  at  outlet,  and  periodicity  or 
connection  conditions  between  the  different  channels  of  the  blade  cascade.  The  time  step,  mesh  and 
turbulence  model  are  chosen  to  put  the  attention  on  the  low  frequency  fluctuations  of  the  attached  cavity, 
more  than  on  the  local  unsteadiness  in  the  cavitation  sheet  wake  (cloud  shedding). 


Cut  at  constant  radius 


Figure  8.1 :  Cascade  and  Mesh  Topology. 

A  non  cavitating  mesh  convergence  was  performed  on  the  following  grids:  mO:  329*61,  ml:  625*121, 
m2:  625*101,  m3:  625*81,  m4:  625*61,  m5:  625*41,  m6:  585*81  and  m7:  841*201.  The  first  number  is 
the  axial  number  of  cells  and  the  second  number  is  the  number  of  cells  in  the  azimuth  direction.  We 
concluded  that  mesh  convergence  was  obtained  with  m4  grid. 

8.2  Steady  Cavitating  Results 

First,  we  present  the  computation  of  head  drop  at  nominal  flow  rate. 
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Figure  8.2:  Head  Drop  Curve. 

The  non  dimensional  head  is  defined  as  the  static  pressure  rise  between  inlet  and  outlet  sections  divided  by 
the  non  cavitating  value.  The  last  converged  point  corresponds  to  a  head  drop  of  17  %. 
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Figure  8.3:  Evolution  of  the  Cavity. 


As  long  as  the  cavity  length  remains  lower  than  half  of  the  distance  between  the  leading  edge  and  the 
throat,  the  head  drop  remains  lower  than  1  %.  Moreover,  the  brutal  head  drop  appears  only  when  the 
cavity  reaches  the  throat.  This  result  will  be  analyzed  thereafter.  Lastly,  it  will  be  observed  that  the 
appearance  of  cavitation  on  pressure  side  takes  place  at  the  same  instant,  which  is  coherent  with  the  results 
that  we  obtained  on  other  geometries  (SHF  impeller  for  example).  Hydraulic  and  mechanical  power  drops 
are  almost  simultaneous  as  it  was  already  observed  on  SHF  impeller  and  on  the  two  inducer  geometries. 
The  loss  of  efficiency  reaches  5,5  %  for  these  simulations.  As  for  the  3D  geometries,  we  examine  the 
evolution  of  blade  load  for  various  a  values: 


Figure  8.4:  Evolution  of  the  Blade  Load. 
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The  cavitation  development  corresponds  to  the  evolution  from  blue  to  red.  The  performance  drop  of  the 
cascade  is  directly  related  to  the  blade  loading:  the  appearance  of  cavitation  on  the  pressure  side  leads  to  a 
torque  decrease  and  then  to  the  performance  drop.  With  regard  to  the  effects  of  compensation,  they 
intervene  here  in  a  rather  significant  way  downstream  from  the  cavity. 

8.3  Unsteady  Cavitating  Results 

In  this  part,  we  present  our  simulations  of  cavitation  instabilities. 

8.3.1  Calculation  Parameters 

The  mesh  that  has  been  used  for  unsteady  cavitating  calculations  is  a  rather  coarse  grid  (215*31)  in  order 
to  minimize  the  computing  time  ;  the  error  induced  is  about  10  %  on  absolute  non  cavitating  pressure  rise 
but  it  is  negligible  regarding  the  prediction  of  head  drop  sigma  value.  Four  channels  have  been  meshed 
since  we  intend  to  simulate  dissymmetrical  phenomena: 


Figure  8.5:  Mesh  Topology  for  Unsteady  Calculations. 


8.3.2  Instability  Configurations 

Several  four-channels  computations  are  performed  at  nominal  flow  rate,  for  different  outlet  cavitation 
number  a  varying  from  low  cavitating  conditions  down  to  the  final  performance  drop.  The  corresponding 
head  drop  chart  is  drawn  in  Figure  5  at  nominal  flow  rate. 
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Figure  8.6:  Calculated  Points  of  Operation. 
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Stable  configurations  correspond  to  cavitation  numbers  lower  than  0.65  or  higher  than  0.8  (square  dots  in 
the  figure).  For  a  >  0.8,  sheets  of  cavitation  are  small  and  identical  in  the  4  channels.  For  these  flow 
conditions,  the  performance  of  the  runner  is  only  slightly  affected  by  the  presence  of  vapour.  For  a  <  0.65, 
the  cavitation  sheets  are  much  more  developed,  and  they  are  responsible  for  the  important  performance 
drop  observed. 

Calculations  are  done  in  relative  referential:  sub  or  super  synchronous  configurations  are  determined  by 
the  sheets  visualization.  For  example,  if  the  sheet  configuration  rotates  in  the  same  direction  as  the 
inducer,  the  characteristic  frequency  of  the  instability  is  added  to  the  inducer  rotation  frequency  and  a 
super  synchronous  configuration  is  established.  In  another  way,  sub  synchronous  regimes  are 
characterized  by  cavitation  patterns  that  rotate  in  the  opposite  direction  as  the  inducer. 

83.2.1  Super  Synchronous  Configuration 

The  super  synchronous  configuration  appeared  for  a  aupstream  value  of  0,75. 


Figure  8.7:  Evolution  of  the  Cavitation  Sheets. 


The  super  synchronous  configuration  is  characterized  by  the  rotation  of  a  small  cavity  in  the  rotating 
direction  of  the  inducer. 

83.2.2  Sub-Synchronous  Configuration 

The  sub-synchronous  configuration  appeared  for  a  aupstream  value  of  of  0,65;  when  a  cavity  was  big 
enough  to  block  a  blade  to  blade  channel: 
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Figure  8.8:  Evolution  of  the  Cavitation  Sheets. 


The  sub  synchronous  instability  is  characterized  by  the  rotation  of  a  big  cavity  in  reverse  direction  from 
the  rotation  of  the  inducer. 

8.3.23  Alternate  Configuration 

An  alternate  configuration  was  obtained  for  ad0wnstream  =  0,75: 


Figure  8.9:  Alternate  Configuration. 

This  configuration  appeared  to  be  rather  unstable  as  we  observed  switches  from  alternate  to  super 
synchronous  configuration.  Another  alternate  configuration  was  found  for  Gd0wnstream  =  0,70  but  this  last 
was  unsteady  in  the  relative  frame. 

8.3.3  Flow  Analysis 

We  devote  this  last  paragraph  to  the  proposal  of  mechanisms  allowing  to  explain  the  development  of 
instabilities.  We  first  study  the  sub  synchronous  configuration  and  then  the  super  synchronous 
configuration. 

8.3.3. 1  Sub -Synchronous  Configuration 

The  appearance  of  a  sub  synchronous  configuration  is  associated  to  the  blockage  of  a  blade  to  blade 
channel,  which  is  quite  close  to  the  stall  phenomenon  observed  at  the  inlet  of  compressors.  Thus,  it  seems 
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to  us  that  the  mechanism  explaining  the  sub  synchronous  instability  may  be  based  on  the  modification  of 
the  incidence  at  leading  edge  caused  by  the  fluctuations  of  mass  flow  in  the  channels: 


Figure  8.10:  Propagation  of  the  Large  Cavity  (Sub  Synchronous  Configuration). 

The  obstruction  of  a  channel  by  a  large  cavity  causes  an  increase  in  the  incidence  on  the  upper  blade, 
which  deepens  the  low  pressure  area  on  this  blade;  the  direct  consequence  is  an  increase  of  the  cavity 
located  on  the  upper  blade.  In  addition,  the  mass  flow  conservation  imposes  an  incidence  decrease  on  the 
blade  below  the  blocked  channel,  which  tends  to  decrease  the  size  of  its  cavity. 

83.3.2  Super  Synchronous  Configuration 

In  order  to  understand  the  propagation  mechanism  of  the  super  synchronous  configuration,  we  use  the 
identification  of  the  blades  presented  on  figure  8.10. 

The  super  synchronous  configuration  is  more  complex  than  the  sub  synchronous  configuration.  From  our 
calculations,  we  propose  an  original  mechanism  based  on  a  coupling  through  the  pressure  field. 
Consequently,  we  observe  the  evolution  of  each  blade  load. 
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Figure  8.1 1 :  Evolution  of  the  Blade  2  Load. 
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The  load  of  blade  2  is  mainly  modified  on  pressure  side.  Thus,  it  is  noted  that  the  pressure  decreases 
gradually  on  first  half  of  the  blade.  This  phenomenon  occurs  whereas  the  cavity  located  in  the  channel 
located  below  is  growing.  This  result  can  be  explained  by  the  development  of  the  low  pressure  area 
associated  with  the  cavity  as  it  was  observed  in  steady  calculations. 

The  load  of  blade  3  is  mainly  modified  on  suction  side  due  to  containment  effects  (decrease  of  blade  4 
cavity  length). 


-3500 

-3580 

3660 

-3740 

-3820 

-3900 


Corde  adimensionnelle 


Figure  8.12:  Evolution  of  the  Blade  3  Load. 


The  blades  concerned  with  a  fluctuation  of  cavity  length  present  the  most  important  loading  variations  as 
pressure  profiles  on  both  sides  are  modified. 


Figure  8.13:  Evolution  of  the  Blade  1  Load. 
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The  growth  of  the  cavity  on  suction  side  suppresses  the  brutal  pressure  rise  located  in  the  closure  area  of 
the  cavitation  sheet.  On  the  pressure  side,  a  peak  appears  near  the  leading  edge  which  is  the  signature  of  an 
incidence  modification.  We  think  that  it  may  cause  the  decrease  of  the  cavity  located  on  blade  4. 
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Figure  8.14:  Evolution  of  the  Blade  4  Load. 


From  the  analyses  of  the  whole  unsteady  results,  we  can  propose  the  following  scenario  to  explain  the 
cavitation  super  synchronous  instability: 

•  The  growth  of  the  cavitation  sheet  at  the  blade  1  suction  side  leads  to  the  increase  of  the  angle  of 
attack  of  this  blade. 

•  A  consequence  of  this  phenomenon  is  the  rise  of  the  pressure  at  the  blade  1  pressure  side,  near  the 
leading  edge. 

•  The  pressure  field  in  the  channel  1,  located  below  the  blade  1,  is  modified  by  this  pressure  rise, 
which  leads  to  the  decrease  of  the  cavitation  sheet  attached  to  the  blade  4. 

•  The  presence  of  the  cavitation  sheet  in  the  channel  4,  located  below  the  blade  4,  will  change  the 
blade  4  angle  of  attack,  inducing  the  progressive  development  of  the  attached  sheet  cavitation,  and 
the  continuation  of  the  instability  propagation. 
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Figure  8.15:  Mechanism  of  the  Super  Synchronous  Configuration. 
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